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ABSTRACT 

We present global three-dimensional self-gravitating smoothed particle hydrodynamics 
(SPH) simulations of an isothermal gaseous disc interacting with an embedded planet. 
Discs of varying stability are simulated with planets ranging from 10 Earth masses 
to 2 Jupiter masses. The SPH technique provides the large dynamic range needed to 
accurately capture the large-scale behaviour of the disc and the small-scale interaction 
of the planet with surrounding material. Most runs used 10 5 gas particles, giving us 
the spatial resolution required to observe the formation of planets. We find four regions 
in parameter space: low-mass planets undergo Type I migration; higher-mass planets 
can form a gap; the gravitational instability mode of planet formation in marginally 
stable discs can be triggered by embedded planets; discs that are completely unstable 
can fragment to form many planets. The disc stability is the most important factor 
in determining which interaction a system will exhibit. For the stable disc cases, our 
migration and accretion time-scales are shorter and scale differently from previously 
suggested. 

Key words: accretion, accretion discs - hydrodynamics - methods: Af-body simu- 
lations - planets and satellites: formation - planetary systems: formation - planetary 
systems: protoplanetary discs 



1 INTRODUCTION 

Planet formation is believed to occur in discs of gas and 
dust around young stars. Observations of circumstellar discs 
in multiple wavelengths reveal the initial conditions for 
planet formati on jMcCaughrean. Stapelfeldt fc Closelboori 
IWilner fc Lavl |200(J) . Detection of extrasolar planets has 
provided a diverse set of final states, including a sta- 
tistically significant number of planets with masses com- 
parable to that of Jupiter on orbits very close to their 
parent star, so-called 'Hot Jupiters' (for a review see 
iMarcv fc ButleJ l200fj . and exoplanets.org for an up-to- 
date catalogue). Formulating an acceptable theory of giant 
planet formation in such an environm ent has proved difficult 
iWuchterl. Guillot fc Lissauerl200(5l . and references therein) . 
Instead planet migration has become the standard theory to 
explain these interesting companions. 

The interaction between a gaseous disc and an embed- 
ded protoplanet is a complex, highly non-linear process. Sev- 
eral phenomena have been classified, including the accretion 
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of planetesimals, the accretion of gas, the exchange of an- 
gular momentum via tidal forces and planet formation via 
gravitational instability. How these behaviours combine is 
determined by the initial conditions, for example the mass 
and temperature of the disc and the location of an embed- 
ded planet. We can sort all possible interactions into a few 
categories: linear migration via tidal interaction, gap forma- 
tion and viscous migration, and planet formation via grav- 
itational instability. All these interactions have previously 
been investigated separately. Here we examine how the in- 
teractions are related and how the initial conditions deter- 
mine which behaviour will be exhibited. 

The dynamics of low-mass planets in gaseous discs have 
been explored both theoretically and numerically. The ex- 
change of orbital angular momentum between disc and em- 
bedded satellite via spiral density waves excited at Lind- 
blad re sonances wa s discussed in iGoldreich fc Tremainj 
(1980). IWardl (Il997f) gave arguments suggesting that this 
will almost always lead to inward migration of the em- 
bedded planet, and derived an expression for the rate 
of this migration in the case where a planet does not 
open a gap in the disc (Type I migration). More recently, 
ITanaka. Takeuchi fc Wardl {2002) gave an updated formula 
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for the total torque on a planet, taking into account coro- 
tation resonances and mentioning three-dimensional (3D) 
effects. Previous numerical simulations have tested some of 
the theore tical assumptions, usually finding genera l agree- 
ment (e.g. iNelson et al"ll2000l : ID'Aneelo et alj|2002h . 

In Type II migration the planet exerts a tidal torque 
greater than the disc viscous force, opening a gap. The 
migration time-scale is slowed, and accretion on to the 
planet is expected to slow or stop. Criteria for gap opening 
have been derived jTakeuchi. Mivama fc Lmlll996t iRafikovl 

2002) and e xamined num e rically using two-dimensiona l 
(2D) grids dBrvden et alJ 1199a ID'Aneelo et all 12002ft . 
The accretion process has been simulated IIKlevI Il999t 
iLubow. Seibert, fc Artvmowiczlllflflflfl finding mass-doubling 
times for Jupiter-mass planets in light discs to be of the or- 
der of 10 5 yr. For a r e view of Type I and II migration, see 
iThommes fc Lissauerl <|2002h . 

The standard picture of migration is situated in a calm, 
stable disc. What if a single planet forms in a ma r ginall y 
stable disc? First hinted at bv lArmitaee fc Hansen! il999t) . 
the spiral density waves excited by a compact object (e.g. a 
rocky core formed via accretion of planetesimals) can lead to 
clump formation that would not otherwise occur. The fact 
that an embedded planet affects the disc density and flow 
locally has been known for some time; here we find that it 
can also have far-reaching (spatially and temporally) global 
effects. 

Regardless of the embedded planets, some discs are 
themselves gravitationally unsta ble. This scenario of planet 
formation has been explored by iBossI lll997t l200lfl and re- 
cently bv lMaver et alJ i2002h . In unstable discs many plan- 
ets of Jupiter-mass or greater can be formed very quickly 
(hundreds of years). The newly formed proto-giant planets 
will violently tear through the disc, and their interaction as 
point masses will dominate the subsequent evolution of the 
system. This scenario is attractive because it easily produces 
planets with properties similar to observed extrasolar plan- 
ets. Furthermore, because the time-scale of formation is so 
small, any system even marginally unstable is likely to form 
planets during the much longer disc lifetime. 

Much of the previous numerical work in this area has 
been limited to two-dimensional approximations, using grids 
with a relatively small spatial range, and ignoring the self- 
gravity of the disc. Some t hree-dimensional simulations have 
been performed recently feate et al.l 120031 ID'Aneelo et alJ 

2003) , although these have still ignored self-gravity. The 
smoothed particle hydrodynamics (SPH) formalism was de- 
signed to be spatially adaptive. Our implementation uses 
variable time-steps, giving us temporal adaptivity. This al- 
lows us to explore the interaction with a much greater dy- 
namic range. Our simulations are fully three-dimensional, 
have free boundary conditions, include the self-gravity of 
the disc and allow the planet to migrate and accrete freely. 

Armed with these useful tools, we have simulated plan- 
ets embedded in gaseous discs, for initial planet masses 
ranging from 10 Earth-masses to 2 Jupiter-masses, and disc 
masses ranging from 0.01 to 0.25 Mq. Previous investiga- 
tions have focused on just one of the interactions described 
above, usually looking at lower-mass discs. Here we attempt 
to relate the interactions, and see how the choice of initial 
conditions determines the classification of the outcome. We 
draw an outline of interesting regions of parameter space, 



and present a few early results from each distinguishable 
region. We also discuss where our simulations differ from 
previous efforts and suggest that this is due to the more re- 
alistic system we are able to model. In future papers we will 
quantitatively address the specific phenomena we can simu- 
late, such as gap opening criteria, gas accretion through the 
gap, Type I and II migration time-scales and the likelihood 
of triggered planet formation. 

The plan of the paper is as follows. In Section we 
describe the initial conditions of our models. The numeri- 
cal technique is briefly discussed in Section [3] and our main 
findings are detailed in Section 3] Details of the results are 
discussed in Section |^| and we conclude with Section [S] 



2 INITIAL CONDITIONS 

We use power-law distributions for density and tempera- 
ture profiles in our simulations, and the disc stability crite- 
rion and mass ratio are dimensionless, so our results should 
scale well. None the less, it is convenient to use conventional 
units when discussing our results. The central star was of 
1 Mq, and its potential was softened with a cubic spline 
on a length-scale of 0.4 AU. All simulations were performed 
in heliocentric coordinates, taking into account the back- 
reaction of the disc and planet on the star. Initially, the disc 
extends from 1 to 25 AU and is free to expand. 

Temperature profile 

Observations of protoplanetary discs do not put tight con- 
straints on the temperature profile. Therefore, we adopt a 
power law for the radial temperature profile and choose 
an index which is equivalent to a constant aspect ratio 
H/r (where H is the scaleheight of the disc at a partic- 
ular radius). Most of our simulations used H/r — 0.05 
(which gives T(5.2 AU) = 102 K). This profile is 
similar to that suggested by the sola r nebula model of 
lHavashi. Nakazawa fc Nakaeawal (I1985T) and is common in 
other simulations of these phenomena. The temperature is 
capped at 5000 K within 0.1 AU of the star and limited 
below by 17 K outside of 30 AU. 

Density profile 

The minimum-mass solar nebula model suggests a power law 
for the vertically averaged surface density of S(r) oc r~ 3 ^ 2 . 
The vertical structure of the disc is determined by the bal- 
ance between the vertical component of the gravity of the 
central star and the vertical pressure support of the disc. If 
we assume no vertical variation in the temperature, the ver- 
tical density structure can be solved, yielding a scaleheight 
expressed as a function of the radius and the temperature 
at that radius. An expression for the density distribution is 
then 

p(r,6,z) = E(r)p vcrt (r,z) = (E 7-- 3/2 ) (p e~ z2/H ^ 

In a gravitationally stable disc we found this profile to 
be stable over many dynamical times, so we chose this 
initial distribution for all of our simulations. Our stan- 
dard disc, 0.1 Mq extending from 1 to 25 AU, thus has 
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£(5.2 AU) = 150 g cm~ 2 . For the gas particles of the 
disc, the initial positions are randomly drawn from the de- 
sired distribution. Initial velocities are chosen so that par- 
ticles are on circular Keplerian orbits, modified slightly to 
account for the self-gravity of the disc and gas forces. In 
stable discs without embedded planets, the surface density 
does not change markedly with time except at the inner 
and outer edges. The outer edge expands slightly, due to 
the initial sharp drop in density. At the inner edge, parti- 
cles migrate inward due to viscous forces, participating in a 
complex balance between the softened potential of the star, 
the smoothed gas pressure forces and Keplerian shear. These 
particles are much closer to the central star than the inner 
Lindblad resonance of the planet, so we do not believe they 
play a large role in determining the evolution of the planet. 
We have not investigated schemes for halting planet migra- 
tion, which might involve removing the inner reaches of the 
disc. 



be several times more massive than the gas particles. There- 
fore, future investigations using higher resolution should be 
able to push the initial planet mass lower. 

2.2 Equation of state 

We use an isothermal equation of state with a spatially fixed 
temperature profile. In doing so we are assuming that the 
disc temperature is controlled solely by the distance from 
the parent star, and that the energy of viscous heating is 
radiated away much faster than the dynamical time-scale of 
the gas. We do not account for shielding, so we miss the 
detailed thermal structure of the clumps that form. Fur- 
thermore, clumps on ell iptica l orbit s will heat up and cool 
down along their orbits. iBossI i2002f) has found that using a 
more realistic equation of state suppresses clump formation 
in marginally stable discs. 



The embedded planet 

The initial mass of the planet particle was a parameter. We 
used planet masses ranging from 10 M® to 2 Mj. The planet 
particle interacts via gravity only, and no gas drag is ap- 
plied. Planets significantly more massive than a gas particle 
can capture gas particles, adding to the effective mass of 
the planet. These captured gas particles are not removed or 
treated specially; they remain as part of the envelope of the 
developing planet. The potential of the planet was spline- 
softened on a length-scale of one-fifth the initial Roche ra- 
dius. 



2.3 Disc stability 

The stability of the disc is characterized by the Toomre Q 
criterion, defined as 



where c s is the sound speed, k is the epicyclic frequency, and 
£ is the surface density of the disc. This criterion measures 
the susceptibility to a local instability. Global spiral waves 
can cause the local density to increase, leading to a local 
collapse. The power laws we have chosen for the temperature 
and density make Q tx r _1,/2 . Hence we expect instabilities, 
if produced, to occur in the outer regions of the disc. 



2.1 An initial gap 

Sufficiently large planets can open a gap in the disc by ex- 
erting a tidal torque stronger than the local viscous angular 
momentum transfer. The balance between these forces de- 
termines the shape and width of the gap formed around such 
a planet, which can be calculated analytically by averaging 
over many orbits and assuming how and where the den- 
sity waves are dampened. The formation of these gaps has 
been studied numerically before, primarily in li ght discs with 
heavy planets without allowing migration fervden et alJ 
1999), finding reasonable agreement with analytic models. 
For low-mass discs, the Type I migration time-scale may 
be sufficiently long for the standard explanation of gap for- 
mation to work. We found that planet migration and gap 
formation could not be decoupled. Therefore, we chose not 
to impose any initial gap, instead the formation of the gap 
is simulated in systems that produce one. 

A Jupiter-mass planet embedded in an unperturbed 
disc without a gap is somewhat unphysical: how did the 
planet become so massive without affecting the disc? As 
seen in Fig. |2l below, the transition from Type I to Type 
II migration is very sudden, taking only a few orbital times. 
Therefore, the sudden appearance of a large planet in a calm 
disc is not entirely unreasonable, and we can imagine the fol- 
lowing scenario: the planet forms, at a lower mass, further 
out in the disc, and undergoes Type I migration and ac- 
cretion to obtain the location and mass we start with. To 
accurately model the capture of gas, the planet particle must 



3 NUMERICAL REALIZATION 

To model the disc-planet system we used Gasoline, a 
tree-based partic le code originally developed for co smolog- 
ical simulations (IStadell l200ll: IWadslev et alJl2004 . It im- 
plements smoothed particle hydrodynamics, a Lagrangian 
treatment of hydrodynamics for gas. It h as been used pre- 
viousl y on solar system-scale p roblems bv lRichardson et alJ 
fcOOd) and iMaver et alJ fcOOSf) . The code was modified to 
apply a fixed temperature profile to the particles. The equa- 
tions of motion include the gravity and hydrodynamics of 
the gas disc. We do not include effects due to magnetic fields. 

3.1 Boundary conditions 

SPH is a Lagrangian technique, so we do not need to impose 
any boundaries. All the particles are allowed to move where 
the forces dictate. Particles are not deleted or merged when 
they approach the central star or the embedded planet. 

3.2 Simulation parameters 

The tree-cell interaction opening angle we used was 
9 — 0.7. Using a multi-stepping technique, the particle 
timesteps are picked according to the local density (with 
r\ = 0.1), and to satisfy the Courant condition (with 
7/c = 0.4). The largest possible timestep was set to be 
2.5 yr, which is 1/50 the orbital time at the outer edge of 
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the disc. The mean molecular weight of the gas is 2 AMU, 
corresponding to molecular hydrogen. 



3.3 Viscosity 

The SPH formalism treats viscosity differently from the 
Shakura & Sunyaev a prescription popular in finite differ- 
ence schemes. No 'real' viscosity is added to the equations 
of motion. Because so little is known concerning the sources 
of real viscosity, we choose to err on the side of caution by 
introducing as little viscosity as possible into the system, 
in an attempt to model the fluid with as high a Reynolds 
number as possible. However, some artificial viscosity is nec- 
essary to stabilize the integration of the equations of motion. 
Gasoline inclu des the standard M onaghan treatment of arti- 
ficial viscosity (lMonaghanlll992l) . including a Balsar a switch 
that reduces viscosity in regions of strong shear fealsaral 
Il995l) . Using the SPH formulation of the equations of mo- 
tion, we can estimat e the kinematic viscosity v = a c s h / 8 
jNelson et alJll99d) associated with the artificial viscosity. 
This estimate is reasonable when there is little or no shock- 
ing (the presence of shocking adds artificial viscosity lo- 
cally). There is strong shear in the disc, so the Balsara switch 
reduces this by a factor of approximately 5, yielding a vis- 
cous time-scale r = r 2 / v > 10 4 yr everywhere for a run 
with 10 5 particles in a stable disc. Since this is at least 5-10 
times longer than our simulations, we believe that artificial 
viscosity will not seriously affect the dynamics of the disc. 
To confirm this, we performed two additional simulations 
where we halved and doubled the standard value of a = 1. 
The change in migration and accretion rates was approxi- 
mately 5 per cent, which is less than the difference due to 
the random seed used to determine the initial particle loca- 
tions. This indicates that artificial viscosity only affects the 
simulation as a numerical stabilizer, and does not introduce 
spurious physical effects. 



4 RESULTS 

Looking first for planet migration, we start with planets em- 
bedded in light, stable discs. Increasing the disc and planet 
masses increases the rate of Type I migration, and brings 
us to gap formation and Type II migration. Continuing to 
marginally stable discs, we find that the density waves ex- 
cited by an already formed planet can trigger gravitational 
instability in discs that would otherwise remain stable. Fi- 
nally, we push the disc mass high enough so that no external 
trigger is needed to drive the gravitational collapse mode of 
planet formation. 

Some of the factors that determine which category a sys- 
tem will end up in include the disc density profile and total 
mass, disc temperature, initial planet mass and orbital ra- 
dius. We made choices for the density and temperature pro- 
files, and varied the total disc mass, the initial planet mass 
and the orbital radius. We ran a simulation for each varia- 
tion of the parameters for a given amount of time (326 yr), 
and classified the final result. We found the disc stability to 
be the most influential factor in determining the behaviour 
of a disc-planet system. The qualitative results of our suite 
of simulations can be seen in Fig. We have arranged our 
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Figure 1. A parameter space map of disc-planet interactions. Q 
is the initial value of the Toomre criterion at the initial location 
of the planet. M p is the mass of the planet, measured in Jupiter- 
masses. The points show the parameter values for the simulations 
we ran. See the discussion in the text. Since the totally unstable 
discs do not require planets, we cannot technically place them on 
this plot. However, if the disc is unstable, an additional planet will 
certainly not prevent instabilities from growing, so it is reasonable 
to label a region of this plot as unstable. 



results in a map of parameter space, illustrating the relative 
positions of the four interactions. 

In stable discs (the upper region) , planets of all masses 
exhibit Type I migration. In more massive discs, planets can 
open a gap, with a bias toward more massive planets. A pre- 
cise measurement of the interface between these two regions 
can show us the criteria for gap opening. To differentiate 
between the regions labelled 'Type I' and 'Gap Formation' 
we look for the halting of the migration of the embedded 
planet within the time we simulated. See Sec. 15.41 for our 
definition of a gap and the reasoning behind it. 

Some planets can trigger the formation of more plan- 
ets via gravitational instability. This requires a marginally 
stable disc and a massive planet. The precise shape we have 
given this region is a guess based on our experience. In the 
bottom region, the disc is gravitationally unstable, even in 
the absence of embedded planets. 

Note that the location of the boundaries between re- 
gions in this map are dependent on the time period of in- 
terest. We are deliberately simulating for only a few hun- 
dred years, to match the time-scale of triggered and wholly 
unstable planet formation. As such, many of the systems 
we classify as Type I do satisfy the standard thermal and 
viscous criteria for gap formation, and would form deeper 
depressions in the density profile and eventually shift from 
Type I to Type II migration if we evolved them several times 
longer. In addition, classification over such a short period of 
time is important, as we measure very high migration and 
accretion rates. 



Simulations of gaseous disc- embedded planet interaction 5 



Table 1. Migration time-scales and accretion rates for Type I 
systems 



M p 


M d (Mq) 


-a (AU/yr) 


M (M Q /yr) 


2 Mj 


0.01 


1.2 x 10" 3 


5.3 x 10" 6 


1.75 Mj 


0.01 


1.2 x 10" 3 


4.7 x 10" 6 


1.5 Mj 


0.01 


1.4 x 10" 3 


4.1 x 10" 6 


1.25 Mj 


0.01 


1.1 x 10~ 3 


3.2 x 10" 6 


1 Mj 


0.01 


1.0 x 10" 3 


2.6 x 10" 6 


0.75 Mj 


0.01 


1.4 x 10" 3 


1.8 x 10" 6 


0.5 Mj 


0.01 


1.1 x 10" 3 


1.2 x 10~ 6 


0.25 Mj 


0.01 


8 x 10~ 4 


2.7 x 10~ 7 


1 Mj 


0.0035 


2.7 x 10" 4 


6.4 x 10" 7 


1 Mj 


0.02 


2.9 x 10" 3 


5.7 x 10" 6 


1 Mj 


0.03 


3.7 x 10~ 3 


8.5 x 10~ 6 


1 Mj 


0.04 


6.0 x 10~ 3 


1.3 x 10~ 5 


1 Mj 


0.05 


7.2 x 10" 3 


1.9 x 10~ 5 


1 Mj 


0.06 


1.0 x 10~ 2 


2.3 x 10" 5 



M p is the planet's initial mass, Afj is the total disc mass, a is the 
inward migration rate, M is the accretion rate on to the planet. 
All these planets started out on circular orbits at 5.2 AU. The 
rates given are equilibrium values over at least the last 200 years 
of the simulation (which typically lasted 326 years). For the lighter 
discs the long-term values were reached immediately. We give val- 
ues for particular simulations, so can give no meaningful error 
estimate. However, we completed several runs at different resolu- 
tion, and usually found agreement within 10 per cent. 




35 70 105 140 175 210 245 280 315 
Time (yr) 

Figure 2. Gap formation scenario: The planet mass as a function 
of time in a 0.1 Mq disc, for different starting masses. The solid 
line is for an initial planet mass of 1 Mj. The top dashed line is 
2 Mj and the bottom dashed line is 0.5 Mj. The dotted line is 
for a 1 Mj planet in the 2D approximation. 



We have shown a two-dimensional slice through the 
three-dimensional parameter space we explored. While con- 
venient and, we believe, largely accurate, this does hide some 
of the complexity we found. For example, we looked for trig- 
gered instability by placing planets at different orbital radii 
in different mass discs such that the value of Q was the same 
at the location of each planet. Since Q is dimensionless and 
we used power laws for the density and temperature, we 
might expect the same result in each system. However, the 
system with closer planet and more massive disc did not 
fragment, while the further planet and less massive disc sys- 
tem did. This may be due to edge effects or a non-linear 
dependence on the effectiveness of Q as a predictor of col- 
lapse. 



as the simulation progresses and a gap is formed. We find 
that the accretion rate scales linearly with the planet mass 
and disc surface density. Further, the migration rate scales 
linearly with the disc surface density. However, changing 
the planet mass has almost n o effect on the migra tion 
rate, confirming the findings of iNelson fc Ben3 l)200otl . If 
we com pare our results wit h the analytic expression (eq. 
70 from lTanaka et ail ^20023), we find our migration time- 
scales [tm = a /H) to be approximately 6 times shorter. 
Most previous simulations of Type I migration have kept 
the planet on a fixed circular orbit, and estimated migra- 
tion time-scales by measuring the torque on the planet from 
the gaseous disc. Our migration timescales are actual mea- 
surements of the rate of change of the orbit of the planet. 



4.1 Type I migration 

Linear theory suggests that a seed planet in a light gas disc 
will undergo Type I migration with a rate that scales linearly 
with the planet mass and disc surface density. The planet 
excites density waves in the disc which torque back on the 
planet, leading to a change in the orbital angular momen- 
tum of the planet. As a result of the power- law index we 
chose for the density profile, we always see inward migra- 
tion. During Type I migration the planet accretes gas from 
the disc with the help of the spiral density waves it excites. 
The internal waves take angular momentum from the disc, 
and the external waves give angular momentum back to the 
disc. On both sides, mass on nearby orbits is shepherded on 
to the planet. 

Tabled lists migration and accretion rates for the sys- 
tems we classified as Type I migration. For the lighter discs, 
the migration rate is a constant function of time, whereas 
a more massive disc causes the migration rate to decrease 



4.2 Gap formation 

More massive seed planets exert strong enough tidal torques 
to open a gap in the disc. By clearing the gas from the loca- 
tion of the Lindblad resonances, the torque exerted on the 
planet by the disc vanishes. Previous theoretical work and 
simulations have suggested that formation of a gap drasti- 
cally reduces the accretion rate on to the planet. We have 
not found this to be true. Even when the gap has cleared, 
the planet still perturbs the disc. These perturbations act 
to pull mass off the edges of the gap and on to the planet. 

We found that gap formation was dependent more on 
the disc mass than the planet mass. This agrees with our 
claim that the self-gravity of the disc is important. At the 
start of the simulation the planet rapidly collects all the 
gas particles nearby, forming a gap within approximately 
100 years. This accretion of gas at corotation is much faster 
than the accretion from close spiral arms or from gap edges, 
as described above. Since the migration time-scale is so short 
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35 70 105 140 175 210 245 280 315 
Time (yr) 

Figure 3. Gap formation scenario: The planet orbital radius as 
a function of time in a 0.1 Mg disc, for different starting masses. 
The solid line is for an initial planet mass of 1 Mj . The top dashed 
line is 2 Mj and the bottom dashed line is 0.5 Mj. The dotted 
line is for a 1 Mj planet in the 2D approximation. 

in the massive discs, there is significant migration during 
the gap formation. After the gap has formed, the migration 
stops or slows considerably. In 0.1 Mq discs, the accretion 
rate slows to a steady value of approximately 10 -5 M0/yr, 
independent of the mass of the planet. Fig.|2]shows the mass 
of a planet and its accreted envelope as a function of time, 
for 0.5, 1, and 2 Jupiter-mass seed planets. By comparing 
with Fig. |21 the semi-major axis of the planet as a func- 
tion of time, we see the accretion slow as the planet clears 
its gap (changes from Type I to Type II migration). The 
planet peels gas off the edges of the gap, as it catches up 
with particles on the outer edge and is caught by particles 
on the inner edge. Our results confirm that the formation of 
'Super- Jupiters' is possible in massive discs. As Type I mi- 
gration ceases, the planet is left in the gap with a non-zero 
eccentricity. We saw eccentricities of 0.01 to 0.04, compara- 
ble to the current eccentricity of Jupiter. 

The time-scale for Type II migration is the viscous time- 
scale of the disc. Since we are actively trying to reduce the 
viscosity of the disc in our simulation, we do not see signifi- 
cant migration after the gap has formed. 



4.3 Planet— triggered instability 

If we reduce the stability of the disc, the spiral arms gen- 
erated by an embedded protoplanet can trigger the gravi- 
tational instability. Such triggered clumps are always exte- 
rior to the seed planet that drives the instability, because 
Q decreases with distance from the star. The spiral arms of 
the planet must present a large enough density contrast to 
drive the local value of Q below 1. Since the unperturbed 
disc stability decreases with distance from the central star, 
for a fixed disc stability, a planet farther out from the cen- 



tral star is more likely to trigger the collapse. If we place 
a Jupiter-mass planet at 12.5 AU, a disc with Q > 1.6 will 
fragment (see Fig. ■ Once the instability is triggered, the 
evolution of the system is very similar to the totally unsta- 
ble case. The clumps form with non-negligible eccentricity, 
and strongly disturb the disc profile. 

4.4 Unstable discs 

If the disc is massive and/or cool enough, gravitational insta- 
bilities will cause it to collapse into bound objects, regard- 
less of any embedded planets. Spiral waves propagate, and 
the density buildup in the arms causes the disc to fragment 
into gravitationally bound objects. Not one, but several of 
these bound objects are formed, and their interaction as 
effective point masses dominates the subsequent evolution. 
The previously smooth disc is ripped apart, as the clumps 
often have sizable eccentricities (see Fig. [5] for a progression 
of snapshots). We observe collisions, ejections, tidal strip- 
ping during close encounters, and violent disruption when 
plunging too close to the parent star. The eccentricities of 
these clumps can be h igh, and change as a result of dynam- 
ical interaction. As in lMaver et al.l i2002f) . we find that the 
disc is unstable when Toomre's Q criterion is less than ap- 
proximately 1.4. The clumps form when a spiral arm with 
sufficient amplitude reaches the outer regions of the disc. 

For this regime we present a simulation that was unsta- 
ble, a 0.2 Mq disc which has Q m i n = 0.8. In this run, a spiral 
arm at approximately 20 AU fragments into four separate 
clumps after 153 years. Five more clumps form over the next 
200 years. At this point, two clumps have developed strong 
leading arms, which fragment into approximately 10 more 
clumps over the next 30 years. Figs.|S]andQshow the orbital 
and mass progression of some of the first clumps formed (the 
simulation forms many more clumps than detailed in the 
plots). We see large eccentricities and eccentricity changes, 
and mergers of several clumps. Once a significant fraction 
of the gas has been caught up in clumps, we expect that 
the furious production of new bound objects will cease, and 
the system will evolve as a set of point masses. After a long 
time (much longer than we have simulated) the system will 
reach dynamical stability, presumably with only a handful 
of objects remaining. 

When a spiral arm fragments, it collapses to a rotating 
spheroid. All our planets formed initially with prograde ro- 
tation, however a later merger resulted in a counter-rotating 
clump. We are not modeling the cooling of the atmosphere 
of the planet accurately, so cannot comment on what the 
final spin of the planets will be. Because our simulations are 
three-dimensional, we can and do observe tilted planets, a 
result of interactions with other planets and vertical asym- 
metry in the spiral density waves. 

Our work on wh olly unstable discs extends the work 
of lMaver et alJ i2002T) to a larger region of parameter space, 
following the system for a longer period of time and spatially 
fixing the temperature profile. The critical value of Q that 
heral ds collapse is quite sensitive to the equation of state 
used (Boss 2002). A more detailed equation of state and 
more physical treatment of radiative transfer are necessary 
to determine how likely planet formation via gravitational 
instability is in a given disc. Note that, though we have 
evolved the system well past the initial collapse, we do not 



Simulations of gaseous disc- embedded planet interaction 7 




Figure 4. Triggered Planet Formation: A sequence of snapshots showing the gas density in a protoplanetary disc. At the planet's 
location Q = 2.2, at the outer edge Q = 1.6. Early, the planet drives spiral density waves in the outer regions of the disc. One of these 
waves increases the local density to the point where the region collapses. The planet continues to drive the spiral arms, leading to further 
clump formation. The process will continue until the disc material is sufficiently depleted. The width of each image represents 55 AU. 



believe we have accurately modeled the detailed structure of 
the planets that form because we are using a simple equation 
of state. However, we have accurately followed the gravita- 
tional interaction of the multiple planets that form. We can 
therefore assert that, in planetary systems that formed via 
gravitational collapse, the disc torque theory of migration 
does not play a role, as the two-body interactions completely 
dominate the orbital evolution of the planets. 



5 DISCUSSION 

Here we describe how we calculated the mass of the planet, 
some tests we performed, further discussion of our migration 
results, and how we define and use the word 'gap' in this 
work. 



5.1 Mass of the planet 

The gas particles that are captured by the planet are not 
removed from the simulation or added to the planet particle. 
They continue to orbit in the disc along with the planet. 
When we give results for accretion rates we have to calculate 
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Figure 5. Unstable Disc Planet Formation: A sequence of snapshots showing the gas density in an unstable disc (Q mln = 0.8). 
Early on, the spiral arms grow. After 150 years one of the spiral arms fragments into a few clumps. As time goes by, other arms fragment, 
in the outer disc. Later, some clumps' leading arms cause fragmentation of the inner disc. The width of each image represents 55 AU. 



a mass of the whole planet, including the initial point mass 
and the captured gas. We use a simple criterion to determine 
whether or not a gas particle counts toward the mass of the 
planet: a particle must be within the Hill sphere of the planet 
and have density greater than the Hill density (the planet 
mass divided by the volume of the Hill sphere) . Since the Hill 
sphere depends on the planet mass, this process is iterated 
until the planet no longer grows, starting from the initial 
planet particle. 



5.2 Resolution issues 

Our standard run for Type I migration is a 0.01 Mq disc with 
a Jupiter-mass planet starting on a circular orbit at 5.2 AU. 
We simulate the system for 326 years (approximately 28 or- 
bits) and calculate accretion and migration rates as a func- 
tion of time. Over this time period the migration and ac- 
cretion rates are constant, and the planet migrates inward 
approximately 0.3 AU and grows in mass to approximately 
1.8 Mj. We ran this model at several resolutions to check for 
convergence. At resolutions of 20000 to 400000 particles the 
accretion rates matched to within 5 per cent and the migra- 
tion rates to better than 1 per cent. Given this agreement, 
we chose 10 5 particles as our standard resolution. 
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Figure 6. Formation in an unstable disc: The orbital radius of 
several clumps as a function of time. The lines begin when a clump 
forms. Clumps that will merge in the time shown have the same 
line style. 
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Figure 7. Formation in an unstable disc: The mass of several 
clumps as a function of time. The lines begin when a clump forms. 
Clumps that will merge in the time shown have the same line 
style. 



The results of two 10 s particle runs can be seen in Fig.|S] 
Each simulation initially has E(r) oc r -3 ^ 2 , aspect ratio 
H/r — 0.05, and an embedded Jupiter-mass planet ini- 
tially on a circular orbit at 5.2 AU. The two density pro- 
files shown are for disc masses of 0.1 and 0.01 Mg, after 
326 years. The disc is three-dimensional, has free boundary 



conditions, and self-gravity is included. The density profiles 
can be compared to the results of prev ious simulation s, for 
examp le fig. 3 of iNelson fc Bend J2003D and fig. 10 of lKlevI 
(1999). Those simulations are both two-dimensional, treat 
viscosity differently, and apply boundary conditions at the 
inner and outer edges of the disc. 



5.3 Migration rate differences 

The Type I migration rates we found were higher than 
most previous numerical and semi-analytical studies. How- 
ever they agree fairl y well with a recent simulation of 
INelson fc Bend i2003t) . which does include the self-gravity 
of the disc. 

We believe our treatment of the vertical dimension and 
the disc self-gravity is correct, and so the discrepancy should 
reveal invalid assumptions made in the linear theory of mi- 
gration. Self-gravity always acts to increase the magnitude of 
density perturbations, so we would expect that the torques 
arising from these perturbations would be correspondingly 
increased, but we have not completed any simulations with- 
out self-gravity, so cannot comment further. To determine 
the effect of including the vertical dimension of the disc, we 
did several runs where all particles were constrained to the 
2 = plane. For light discs, the 2D case consistently gave ac- 
cretion rates lower than the 3D case by at least 20 per cent. 
For more massive discs, the 2D case appears to give a slightly 
higher accretion rate during gap formation but again a lower 
rate after the gap has formed (see Fig. The difference in 
migration rates between 2D and 3D simulations was not 
clear. 



5.4 Defining a gap 

Determining if and when a planet forms a gap in a disc is a 
non-trivial problem. By what percentage must the density 
decrease around a planet to be called a gap, and how far 
must this density depression reach? These criteria necessar- 
ily invoke an arbitrary number. Fig.|H]shows an azimuthally 
and vertically averaged surface density of two simulations 
with different disc masses. In both cases the gas density 
near the planet has been markedly reduced, and these de- 
pressions could be called gaps. The planet in the heavier 
disc has stopped migrating (see Fig. |3J. In the lighter disc, 
however, the planet continues to migrate inward at a rate in- 
dependent of the depth of the density depression. Therefore, 
we choose to use the word 'gap' only when its dynamical ef- 
fects are apparent, i.e. when the planet has stopped migrat- 
ing. This is how we chose between the 'Type I' and 'Gap 
Formation' regions when classifying simulations for Fig. 



6 CONCLUSION 

This work presents the first results from our ongoing investi- 
gation of planet formation and migration. It is a pioneering 
application of global SPH simulations to a field previously 
examined primarily with Eulerian grid-based simulations. 
Using the new technique, we have measured migration and 
accretion time-scales, witnessing Type I migration and gap 
formation. Our results for this region partially qualitatively 
agree with the linear theory of interaction with spiral density 
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Figure 8. Azimuthally and vertically averaged surface density of 
two discs with embedded Jupiter-mass planets, after 326 years. 
The solid line is the density of a 0.1 Mq disc. A deep gap has 
formed within 100 years, and Type I migration has stopped. The 
dotted line is the density (multiplied by 10) of a 0.01 Mq disc. 
Note that a depression has formed about the planet, of less than 
half the initial density. This feature can be called a gap, but the 
planet is still exhibiting Type I migration. The dashed line shows 
the initial density profile. 

waves. We have observed planet formation via disc fragmen- 
tation, including a cascade effect from a single starter planet. 
Using massive discs, we find that it is easy to form large (tens 
of Jupiter masses) planets on a short time-scale (hundreds 
of years). If planets form via gravitational instability, the 
theory of migration via disc torques is not applicable. 

To obtain an overview of gaseous disc-embedded planet 
interactions, we performed a suite of simulations, varying the 
planet mass, the initial location of the planet, and the disc 
stability. We have shown the orientation in parameter space 
of four broad categories of interaction. These interactions 
are: Type I migration; gap formation and Type II migration; 
planet-triggered instability; and unstable discs. The stability 
of the disc is the most effective predictor of the type of 
interaction between a gaseous disc and an embedded planet. 

Planet formation and migration is a non-linear process 
that requires the reasoned application of computational re- 
sources to model. Future investigations will explore the de- 
tails of each interaction we have described. In addition, we 
will refine our knowledge of the boundaries between regions 
of parameter space. Future observations will give us a range 
for Q and a better understanding of the temperature pro- 
file in protoplanetary discs. Combined with the results pre- 
sented here, we can say which migration and formation pro- 
cesses will be relevant. 
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